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O ; ABSTRACT 

, Hot spots residing on the surface of an accretion disc have been considered 

' as a model of short-term variability of active galactic nuclei. In this paper 

we apply the theory of random point processes to model the observed signal 
from an ensemble of randomly generated spots. The influence of general 
relativistic effects near a black hole is taken into account and it is shown 
r-j I that typical features of power spectral density can be reproduced. Con- 

^ |. nection among spots is also discussed in terms of Hawkes' process, which 

produces more power at low frequencies. We derive a semi-analytical way 
to approximate the resulting power-spectral density. 
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> ■ 1 INTRODUCTION 

(N : 

. Radiation from accreting black holes varies on different timescales [1]. In X-rays, the 

' observed Ught-curve, / = /(t), is a complicated noisy curve that can be represented 

, by a broad-band power spectrum Pj. It has been proposed [S] |3] that 'hot spots' 

are a possible contributor to this variabiUty. These spots are supposed to occur on 
the surface of an accretion disc following its irradiation by coronal flares [3 [6l [7] . A 
model light-curve can be constructed as a sum of contributions from many point-like 
sources that are orbiting above an underlying accretion disc. The observed signal is 
modulated by relativistic effects as photons propagate towards a distant observer. 
In order to characterise the light curves we need to introduce some appropriate 
^ • estimator of the source variability. In a mathematical sense, one applies a functional: 

5^ , f ^ S [f], where S [.] is a map from functions defined on R to functions on R'^ 

(fc > 0). The variability estimator can be a single number (for example, the mean 
flux or the 'rms' characteristic), or function of one variable (power spectrum density 
or probability distribution) or of many variables (poly-spectra, rms-flux relation, 
etc). A signal of such a spotted accretion disc should be intrinsically stochastical. 
Hence, the variability estimator S [/] , derived from a piece of the light-curve, is a 
random value, too. 
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Various schemes have been proposed in which spots are mutually interconnected 
in some way [51 [S]- We want to investigate this type of models within a com- 
mon mathematical basis. Other authors have developed different approaches to the 
problem [TUlin]- 



2 SPOT MODELS AND THE ACCRETION DISC VARIABILITY 

2.1 Model assumptions and variables 

Let us have K samples of the observed light-curves from the same source, fj . The 
law of large numbers ensures that 

j^J2^[fj]^^[S[f]], ^^oo, (1) 

where E[.] is the mean value operator. The average value of the functional is formally 
defined 

E[5[/]]= J2 (Probability of /,) X 5 [/,], (2) 

{All possible /j(t)} 

where the sum goes over all possible light-curves generated by this model. We will 
show how to define and parameterise "the space of all possible light-curves" and 
how to perform the averaging when the functional is the power spectrum. 

The general model is constrained only by the following three assumptions about 
the creation and evolution of spots: 

(i) Each spot is described by its time and place of birth (tj , rj and (f>j ) on the 
disc surface. 

(ii) Every new occurrence starts instantaneously; afterwards the emissivity decays 
gradually to zero. The total emitted radiation energy is finite. 

(iii) The intrinsic emissivity can be fully determined by a finite number of pa- 
rameters which form a vector ^j. 

For a simple demonstration of this concept see figured] The disc itself has a passive 
role in our present considerations. We will treat it as a geometrically thin, optically 
thick layer lying in the equatorial plane. 

2.2 Random point processes 

The concept of point processes is a generalisation of well-known random processes 
which were developed as a description of time-dependent random values [12| . Point 
processes are used as statistical description of configurations of some randomly 
distributed points in space M". 
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One way of describing a configuration of points is by their counting measure, 
N{A), which for every A C M" gives a number of points lying in A. One defines 
the intensity measure, 



Ml (A) = E[iV(A)] 



(3) 



Similarly to random processes, the point process can be characterised by its mean 
value and moments. For every A C M", Mi{A) is the mean number of points lying in 
A. The second-order moment measure is defined in the same way on the Cartesian 
product of spaces R" x R": 



M2{A X B)^E [N{A)N{B)] 



(4) 



Let Xipf be one possible configuration of points, i.e. the support of some N{.) For 
the functions f{x) and g{x, y) on R" and R^", respectively, it follows [nifTlj 



E 



{xi}N 



E 



E yi 

{■■i:i}N, {yi}N 



= E 



E 



f{x)N{dx) 



X 



= / /(x)Mi(da:) (5) 



g{x, y)N{dx)N{dy) 
= I g{x, y)M2{dx x dy). 



(6) 



The concept of point process can be further generalised in the following way. We 
add a mark Hi from the mark set K. to each coordinate Xi from {xi}^- Marks carry 
additional information. The resulting point process on the set R" x /C is called 
the 'marked point process' if for every A C R" it fulfills the condition Ng{A) = 
N{A xK.) <oo. 

The random measure Ng{A) represents the ground process of the marked process 
N. When the dynamics of the process is governed only by the ground process and 
marks are mutually independent and identically distributed random values with 
the distribution functions G(d/v), then the process intensity and the second order 
measure fulfill 



Mi{dx X dK) 
M2(da;i x d^i x dx2 x dK2) 



Mig(da;)G(dK), (7) 
Afg2(dxi X dx2)G(dKi)G(dAt2). (8) 



2.3 Relationship between point processes and spots 

Let us assume a surface element orbiting at radius r with constant emissivity / and 
orbital frequency ^{r). This should represent an infinitesimally small spot. For the 
flux measured by an observer at inclination 0^ we find 



f{t)=IF{t, r, 0o). 



(9) 
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The periodical modulation of the signal is determined by relations 

Fiti(t>),r,9o)=F{<P,r,9o), (10) 
m^T^+St{^,r.,e,), (11) 

where r, 9o) is the transfer function describing the total amplification of signal 
emitted from then disc surface element on the coordinates r and (p. The function 
St{(j), r, 9o) is the time delay of the signal (hereafter we will omit 9o in the argument 
of F for simplicity). Now, we consider a process consisting of statistically dependent 
events, 

fit) - 5]/(i - S,, ij)F{t - <5p„ rj), (12) 
j 

where: I{t, ^) = 9{t)g{t, ^) is the profile of a single event; 6j = tj + toj is time 
offset; Spj = 6j + tpj is the phase offset; 6{t) is the Heaviside function; and g{t, ^) 
is non-negative function of k + 1 variables t and ^ = (^^, . . . ,^'^), which is on the 
interval (0, oo) integrable in the variable t for all values of parameters ^ S S. The 
set S is some measurable subset of M*^. For a fixed value of r, F{t, r) is a periodical 
function of t, with the angular frequency ri(r). 

Quantities ^j, tj, rj, tpj and toj are random values. The vector ^j determines 
the duration and shape of each event, tj is time of ignition of the j-th event, and 
toj the corresponding initial time-offset. Parameter tpj determines the initial phase 
of the periodical modulation. Processes of this kind and their power spectra were 
mathematically studied by Bremaud and Massoulie, [TFl I16| . 

Power spectral function of a stationary stochastic process X{t) is 



lim ;5^E \TT[Xmu;)f , (13) 



1 

2^ 

where Tt [ ] is the incomplete Fourier transform. 



rT[X{t)] = j X{t)e~'^'dt. (14) 

-T 

This can be evaluated by using the complete Fourier transform, 

T oo 

X(Oe— *dt= / X(i)x(-T.T>(Oe-'"*di = 25^^^*.F[X(t)](^), (15) 



UJ 



where xa {x) is the characteristic function of set A, which equals 1 for x ^ A and 
for X ^ A. Symbol ★ denotes the convolution operation. 

By applying this transformation on the process (|12p we find 

TT[f{t)]{oj) = ^^^^*^^[/(t~<5,, F{t-6p„ rj)]{u). (16) 



Modeling an accretion disc stochastical variability 5 




10.0 - c) 
8.0 - 
e 6.0 - 



4.0 - 
2.0 - 

0.0 I ' ' ' ' ' ' ' ' ' 

10 20 30 40 50 60 70 80 90 100 

t 

Figure 1. The model light-curve (panel a) is obtained as a sum of elementary events 
(panel b). Profile of each individual event is assumed to be I{t,T) = (t/r)^ exp (—t/r) 9{t), 
as described in the text. Their normalization is identical. The final light-curve is fully 
determined by the form of the individual contributions together with a set of points in t~T 
plane (panel c), which represent pairs of ignition times and temporal constants r of each 
event. 



The Fourier transform of a single event I{t — Sj, ^j)F{t — Spj, rj) is then 

T[I{t ~ 6„ i^)F{t - Sp„ r,)](w) = e-'^-'^J^[I{t, i,)]*nF{t + h,, r^)]- (17) 

Function 7') is periodical in time, and so it can be expanded: 

oo 

F(t, r)= Cfc(r)e^'="«*, (18) 

k— — oo 

where ^{r) is the frequency of _F(t, r). We find 



T[F{tp,r)]{u:)= ^ Cfe(r)e'''^"W*pJ (c. - fcl7(r)) , (19) 

k — — OQ 

oo 

T [I{t, 0] * ^ [F{t + ip, r)] - Cfe(r)e"="('-)*p.F [/(i, 0] - kn{r)) . (20) 



fc— — oo 



The above given formulation of the problem falls perfectly within the mathematical 
framework of point processes. 
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2.4 The case of independent decaying spots (Poisson process) 

Knowing the incomplete Fourier transform of f{t) we can now calculate its squared 
absolute value and perform the averaging over all realizations of the process. Be- 
tween — T and T the process is influenced by all events ignited during the preceding 
interval (—00, T), however (because of fast decay of every single event), this can 
be restricted onto {—{T + C), T), where C is a sufficiently large positive constant. 
Therefore, every realization of the process f{t) on the interval (— T, T) can be de- 
scribed by set of points in (k + 4)-dimensional space {tj, tQj, tpj, rj, ^j), where 
e(-(T + C),T). 

Equation (|12p represents a very general class of random processes. However, in 
all reasonable models of spotted accretion discs the values of initial time delay and 
phase are functions of initial position of each spot (r and (j)), i.e. 



to = St{r, (j)), tp = —— + to- 



(21) 



Fourier transform of the resulting signal can be then simplified, 

- hj, ij) F{t - toj + tp,, r,)]{u;) 

00 

= E c,{r)e^''^T[I{t-5t{r,cj,),^)]{u;-kn{r)). (22) 



k— — c 



Every realization of this process is completely determined by set of points {tj , cjjj , rj , ] 
from some subset of R'^^^ . 

For the sum of K complex numbers Zi it follows 



K 



K 



K 



K 



K 



K K 



i=l ] = ! 



\i=l 



Defining the function s{t,(f>,r,^;u!) as 
s(t, (j), r. 



|;^)^i£!!^*|e--t E cu{r)e^^^:F\I{t~5t,£,)]{u~kn{r)) 

\ k——-oo / 



According to (f23j) we can write 



(23) 



(24) 



I^T[/(t)](c.)| 



Es(ij, if ^) 



= EE'^(*J' '^J' ''j' ^3' '^)^*(^'' ^i' "^i' ii' (25) 

j I 
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E 



\^T[fit)]iu;)\' 



E 



'^J' ^j' ^) "^^J' ^) 

sit, 0, r, e 0', r', a>)m2(i, 0, r,^ t' , 0', /, OdAdA', (26) 



AxA' 



where m2 is density of the second-order moment measure corresponding to the 
random point process of (tj, 4>j, r^, ^j). The set ^ is a Cartesian product of sets, 

A = (-(T + C), T) X (0, 27r) x (w, w) x S. (27) 

Now we can perform the limit (|13p . It can be shown that the result is independent 
on the value of C. In order to obtain an explicit formula for the power spectral 
density we have to specify the form of M2(.). In the simplest case we assume 
events that are mutually independent with uniformly distributed ignition times. 
The process can be described as a marked point process with a Poissonian process 
as the ground process. The intensity and the second-order measure for the ground 
process are: 

Mgi(di) = ndt, (28) 
Mg2idtdt') = [n^ + n5{t-t')]dtdt', (29) 

where n is the mean rate of events. Other parameters are treated as independent 
marks with common distribution G(d0drd^). The second-order measure of the 
process has a form 

M2(did0drd|dt'd0'd^') = [n'^G{d(j)dr d$,)G{d(j)' dr' d^') + nG{d(j)dr d^) 

X 5{t - t')S{(j) - (j)')5{r - r')S{^ - ^')] dtdt'. (30) 

For the power spectrum we obtain this general formula, 

oo oo „ 

S{uj) = 4n^n J2 Ck{r)c*i\r)e'^^~^^^ T [I{t - 5t{r, 0), ^)] {lo - kn{r)) 

k— — oo l—^oo^ 

X J-* [I{t - (5i(r, 0), ^)] {oj - m{r)) G'(d0drd^). (31) 



2.5 Introducing a relationship among spots (Hawkes process) 

The assumption that the spots are mutually statistically independent seems to 
be a reasonable first approximation, however, the actual ignition times and spot 
parameters should probably depend on the history of a real system. As an example 
of such non-Poissonian process, we calculated the power-spectral density (PSD) for 
a model in which the spot ignition times are distributed according to the Hawkes 
[17] process. 
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The Hawkes process consists of two types of events. Firstly, new events are 
generated by Poisson process operating with the intensity A. Secondly, an existing 
event with ignition time ta can give birth to new event at time t according to Poisson 
process with varying intensity ^{t ~ to). So the mean number of events found at 
time t is 

m(<) = A+ ^ ii{U) = \+ I n{t)N{At). (32) 

i,ti<t 

For a stationary process the first moment density is constant. Averaging both sides 
of the previous equation we find, 

oo 

mi — , v ^ j ii{t)dt. (33) 

— oo 

Stationarity of the process implies, that the second-order measure density can de- 
pend only on the difference of its arguments. It can be proven [14) that 

mg2(t, t') = c(t - t') + ra\^ + m^i8(t - t'), (34) 

where the c(i) is an even function. Thus, for the corresponding marked process 
with independent marks we find M2(dt dc/) dr d^ dt' d0' d^'): 



M-2 = 



\2 

+ c{t-t')] G{d(p dr d|)G(d0' dr' d^') 



+ ^^G(d(j)dr dOSU - t')S(<j) - <f)')S(r - r')S(^ - $') 
1 — 1/ 



dtdt'. (35) 



This second-order measure is almost identical to that of the Poissonian process 
(there is only one additional term associated with the function c{t)). The resulting 
PSD is 



S{uj) = 47r2^ ^ ^ y e^^'->^^^Ck{r)c*{r)T [Kt - 5t{r, 0] - ^^^(r)) 

k— — oo l — — oo^ 

X T* [I{t ~ St{r, 0), $,)] {uj - m{r)) G{d(f) dr d|) + in^T [c{t)] (uj) 

OO p 

^ Ck{r) J e~''^J'[l{t-St{r,cP),^')]{u;-kn{r))G{d<t>drd^) 



X 

k— — oo 



oo „ 

^ J e^'^'cUr')^* [I{t - 6t{r', 0'), 0] ^ mr')) G(d0' dr' d^'). (36) 



X 

/ — — oo 



The function c{t) can be calculated from the mean number of secondary events 
fj-(t). Assuming /i(t) = z/a exp(— ai)0(t) we obtain 



Xaiy{l ~ vl2) 



= n ..,2 exp(-a(l - ,y)\t\). (37) 
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Figure 2. Power spectra from the spot model driven by tiie Poisson process (top row) 
and the Hawkes process (bottom row), calculated for spots orbiting and evolving on the 
surface of a thin accretion d ISC i^Tin = 6, Tout ~ 100 gravitational radii). Two values 
of observer's inclination 6o are shown for comparison. The red (thin, noisy) curve is a 
result of direct numerical simulation. Blue (thick, continuous) curves are the analytical 
approximations based on eqs. pip and (|36|l , respectively. We assumed probability density 
function p(r) oc 1/r. The magenta (vertical) lines denote the Keplerian orbital frequency 
r2(r) at the inner and the outer edges of the disc. One can see that the Hawkes' process 
tends to enhance the low-frequency part of PSD and shift the break frequency towards 
lower values, below f2(rout)- 



It is interesting to notice that the above-given formal approach can actually pro- 
vide a useful analytical formula to approximate the power spectrum. Figure [2] shows 
exemplary PSD which were obtained by (i) direct computations of the light-curve 
and the resulting PSD, and by (ii) the semi-analytical approach with Poissonian 
and Hawkes processes. 



3 CONCLUSIONS 

We have studied the properties of power spectral density within the model of ac- 
cretion disc variability driven by orbiting spots. The origin and evolution of spots 
were described in terms of Poissonian and Hawkes' processes. The latter belongs to 
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a category of avalanche models. We developed an analytical approximation of PSD 
and compared it with our numerical results from light-curve simulations. In this 
way we were able to demonstrate the precision of formulae (|3ip and The an- 

alytical approximation evaluates very fast and provides the main trend of the PSD 
shape while avoiding the noisy form of the numerically simulated spectra. Our ap- 
proach allows us to investigate the resulting PSD as a function of the assumed type 
of process, which describes creation of parent spots and the subsequent cascades 
of daughter spots. In particular, we can investigate the predicted PSD slope at 
different frequency ranges and we can locate the break frequency depending on the 
model parameters. 

The resulting PSD can be approximated by a broken power-law. For every sta- 
tionary process the quantities S(0) and S{Lu)du! are finite. Therefore, the func- 
tion S{lu) flattens {S{oj) « uj'^) near to — and it must decrease faster than 1/uj at 
high frequencies. Power-spectra generated by the spot model behave in this way. 
The low- frequency limit is a constant, whereas the-high frequency behaviour de- 
pends mainly on the shape of the spot emission profile, I(t, ^j). In our calculations 
the emissivity was a decaying exponential and the slope was equal to —2 at high 
frequencies. The most interesting part of the spectra in between those two limits is 
influenced by both the emissivity profile and the underlying process. 
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